/********************************************************
 * @file: matrix.h
 * @date: 25/09/2008
 ********************************************************
 * FBElem v1_0 - (c)copyright 2008 GoBEM
 ********************************************************/

#ifndef _MATRIX_H_
#define _MATRIX_H_

// includes
//----------
//#include <util/typedefinitions.h>
#include <iostream>
#include <assert.h>
#include "vector.h"
#include <math.h>
#include <iomanip>

// forward declarations
//----------------------
template <class T> 
class Vector;

// useful typedefs
//-----------------
using namespace std;


/*!
 * CLASS
 * Matrix
 * 
 * This class defines the capabilities of a matrix (a two dimensional array).
 *
 * It has several functions which aid the matrix manipulations and operations,
 * is an efficient way. But one should think in use this as a wrapper class in
 * order to have access to different LIBs.
 *
 * Since it is an template, it allows the use of both double and complex matrices.
 * 
 * @author Andre Pereira (andremaues@gmail.com)
 * @version 1.0.0
 * @see Vector
 */
template <class T>
class Matrix
{

public:

  // Constructors
  Matrix( unsigned int _rows = 2, unsigned int _cols = 2);
  Matrix(const Matrix<T>& m);

  // Destructor
  ~Matrix();

  // Assignment operator
	Matrix<T>& operator = (const T& _a);
	Matrix<T>& operator = (const Matrix<T>& m);

  // Value extraction method
  unsigned int getRows() const;
  unsigned int getCols() const;

   // Subscript operator
  virtual T& operator () (unsigned int row, unsigned int col);
  virtual const T& operator () (unsigned int row, unsigned int col) const;
  // Calculation operators
  virtual Matrix<T> operator + (const Matrix<T>& m) const;
  virtual Matrix<T> operator - (const Matrix<T>& m) const;
  virtual Matrix<T> operator * (const Matrix<T>& m) const;
//  Matrix<T> operator * (const CompressedColumnStorage<T>& m) const;
  virtual Vector<T> operator * (const Vector<T>& v) const;
  virtual Matrix<T> operator * (const T& a) const;

  // Combined assignment - calculation operators
  // Vector<T>& operator += (const Vector<T>& v);

  // Equality operators
  virtual bool operator == (const Matrix<T>& v) const;
  virtual bool operator != (const Matrix<T>& v) const;

  // Miscellaneous -methods
  void zero();
  virtual void zeroRowColumn(unsigned int _k);
  // void unit();

  // Utility methods
	void resize(unsigned int _rows, unsigned int _cols);
	void resize(unsigned int _n);
  void addPartMatrix(	const Matrix<T>& _m,
                      const unsigned int _startRow,
                      const unsigned int _startCol );
  void addPartMatrixWithColVec(	const Vector<T>& _v,
                                const unsigned int _startRow,
                                const unsigned int _startCol );
  void setPartMatrixWithColVec(	const Vector<T>& _v,
                                const unsigned int _startRow,
                                const unsigned int _startCol );
  Vector<T> getColumnVector( const unsigned int _iColNumber ) const;
  Vector<T> cross( const unsigned int _i,
                   const unsigned int _j ) const;
  virtual void print() const;
  const bool equalWithTolerance( const Matrix<T>& _m, const T& _a ) const;
  // det();
  // transpose();
  // inverse();
  virtual Vector<T> solve( Vector<T>& _v );

  // Type of matrices

private:

    T * m_data;
    unsigned int m_rows, m_cols;

};

typedef Matrix<double> DMatrix;
typedef Matrix<int>    IMatrix;


//-----------------------------------------------------
template <class T> inline
Matrix<T>::Matrix( unsigned int _rows, unsigned int _cols )
{
	m_data = new T[ _rows * _cols ];
	assert( m_data != 0 );
	m_rows = _rows;
  m_cols = _cols;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>::Matrix(const Matrix<T>& m)
{
  m_rows = m.m_rows;
  m_cols = m.m_cols;
  unsigned int size = m_rows * m_cols;
	m_data=new T[size];
	assert(m_data!=0);
	for(unsigned int i=0;i<size;++i)
    m_data[i]=m.m_data[i];
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>::~Matrix()
{
  delete [] m_data;
  m_data = 0;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>& Matrix<T>::operator =(const T& _a)
{
    unsigned int size = m_rows * m_cols;
  for( unsigned int i=0; i<size; ++i )
    m_data[i] = _a;
	return *this;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>& Matrix<T>::operator = (const Matrix<T>& m)
{
	if(this != &m)
	{
		delete [] m_data;
    m_rows = m.m_rows;
    m_cols = m.m_cols;
    unsigned int size = m_rows * m_cols;
		m_data= new T[size];
		assert(m_data!=0);
		for(unsigned int i=0;i<size;++i)
       m_data[i]=m.m_data[i];
	}
	return *this;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
unsigned int
Matrix<T>::getRows() const
{
  return m_rows;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
unsigned int
Matrix<T>::getCols() const
{
  return m_cols;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
T&
Matrix<T>::operator () (unsigned int _row, unsigned int _col)
{
	assert( _row <= m_rows && _col <= m_cols );
  return m_data[_row*m_cols+_col];
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
const T&
Matrix<T>::operator () (unsigned int _row, unsigned int _col) const
{
	assert( _row <= m_rows && _col <= m_cols );
  return m_data[_row*m_cols+_col];
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>
Matrix<T>::operator + ( const Matrix<T>& m ) const
{
	assert( m.m_rows == m_rows && m.m_cols == m_cols );
  Matrix<T> sum( m_rows, m_cols );
  unsigned int i, j;
  for( i=0; i<m_rows; i++ )
	  for( j=0; j<m_cols; j++ )
		  sum.m_data[i*m_cols+j] = m_data[i*m_cols+j] + m.m_data[i*m_cols+j];
  return sum;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>
Matrix<T>::operator - ( const Matrix<T>& m ) const
{
	assert( m.m_rows == m_rows && m.m_cols == m_cols );
  Matrix<T> sub( m_rows, m_cols );
  unsigned int i, j;
  for( i=0; i<m_rows; i++ )
	  for( j=0; j<m_cols; j++ )
		  sub.m_data[i*m_cols+j] = m_data[i*m_cols+j] - m.m_data[i*m_cols+j];
  return sub;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>
Matrix<T>::operator * ( const Matrix<T>& m ) const
{
	assert( m_cols == m.m_rows );
  unsigned int i, j, k;
  unsigned int cols = m.m_cols;
  Matrix<T> mult( m_rows, cols );
  for( i=0; i<m_rows; i++ )
	  for( j=0; j<cols; j++ )
    {
      mult.m_data[i*cols+j] = T(0);
      for( k=0; k < m_cols; k++)
		    mult.m_data[i*cols+j] += m_data[i*m_cols+k] * m.m_data[k*cols+j];
    }
  return mult;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T>
Matrix<T>::operator * ( const Vector<T>& v ) const
{
	assert( m_cols == v.getLength() );
  Vector<T> mult( m_rows );
  unsigned int i, j;
  for( i=0; i<m_rows; i++ )
  {
    mult(i) = T(0);
	  for( j=0; j<m_cols; j++ )
      mult(i) += m_data[i*m_cols+j] * v(j);
  }
  return mult;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>
Matrix<T>::operator * ( const T& a ) const
{
  Matrix<T> result(m_rows,m_cols);
  unsigned int size = m_rows * m_cols;
  for( unsigned int i=0;i<size;i++ )
    result.m_data[i] = m_data[i]*a;
  return result;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T>
bool
Matrix<T>::operator == (const Matrix<T>& m) const
{
	if( m.m_rows != m_rows || m.m_cols != m_cols ) return false;
  unsigned int size = m_rows*m_cols;
  for(unsigned int i=0;i<size;i++)
    if(m_data[i] != m.m_data[i])
      return false;
	return true;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T>
bool
Matrix<T>::operator != (const Matrix<T>& m) const
{
    return !this->operator ==(m);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void
Matrix<T>::zero()
{
  unsigned int i, j;
  for( i=0; i<m_rows; i++ )
	  for( j=0; j<m_cols; j++ )
		  m_data[i*m_cols+j] = T(0);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void
Matrix<T>::zeroRowColumn( unsigned int _k )
{
  unsigned int i, j;
  for( i=0; i<m_rows; i++ )
    m_data[i*m_cols+_k] = T(0);
  for( j=0; j<m_cols; j++ )
    m_data[_k*m_cols+j] = T(0);
  m_data[_k*m_cols+_k] = T(1);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void 
Matrix<T>::resize(unsigned int _rows, unsigned int _cols)
{
	delete [] m_data;
  m_rows = _rows;
  m_cols = _cols;
  unsigned int size = m_rows * m_cols;
	m_data = new T[size];
	assert(m_data !=0);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void 
Matrix<T>::resize(unsigned int _n)
{
	delete [] m_data;
  m_rows = _n;
  m_cols = _n;
  unsigned int size = m_rows * m_cols;
	m_data = new T[size];
	assert(m_data !=0);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void
Matrix<T>::addPartMatrix(	const Matrix<T>& _m,
                          const unsigned int _startRow,
                          const unsigned int _startCol )
{
	assert( _startRow < m_rows );
	assert( _startCol < m_cols );
  unsigned int R = _m.m_rows;
  unsigned int C = _m.m_cols;
	assert( ( _startRow + R ) <= m_rows );
	assert( ( _startCol + C ) <= m_cols );
  unsigned int i, j, k;
  for( i=0; i<R; ++i )
  {
    k = i + _startRow;
    for( j=0; j<C; ++j )
	    m_data[k*m_cols+j+_startCol] += _m.m_data[i*C+j];
  }
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void
Matrix<T>::addPartMatrixWithColVec(	const Vector<T>& _v,
                                    const unsigned int _startRow,
                                    const unsigned int _startCol )
{
  assert( _startRow < m_rows );
  assert( _startCol < m_cols );
  unsigned int L = _v.getLength();
  assert( ( _startRow + L ) <= m_rows );
  for( unsigned int z=0; z<L; z++ )
    m_data[(_startRow+z)*m_cols+_startCol] += _v(z);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void
Matrix<T>::setPartMatrixWithColVec(	const Vector<T>& _v,
                                    const unsigned int _startRow,
                                    const unsigned int _startCol )
{
	assert( _startRow < m_rows );
	assert( _startCol < m_cols );
  unsigned int L = _v.getLength();
	assert( ( _startRow + L ) <= m_rows );
  for( unsigned int z=0; z<L; z++ )
	  m_data[(_startRow+z)*m_cols+_startCol] = _v(z);
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T>
Matrix<T>::getColumnVector( const unsigned int _iColNumber ) const
{
	assert( _iColNumber < m_cols );
  Vector<T> col( m_rows );
  for( unsigned int i=0; i<m_rows; i++ )
    col(i) = m_data[i*m_cols+_iColNumber];
  return col;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T> 
Matrix<T>::cross( const unsigned int _i,
                  const unsigned int _j ) const
{
  assert( m_rows <= 3 && m_rows > 1 );
  Vector<T> crossProd(3);
  if( m_rows == 3 )
  {
    crossProd(0) =   m_data[m_cols+_i]*m_data[2*m_cols+_j] - m_data[2*m_cols+_i]*m_data[m_cols+_j];
    crossProd(1) = - m_data[_i]       *m_data[2*m_cols+_j] - m_data[2*m_cols+_i]*m_data[_j];
  }
  else
  {
    crossProd(0) = 0.0;
    crossProd(1) = 0.0;
  }
  crossProd(2) = m_data[_i]*m_data[m_cols+_j] - m_data[m_cols+_i]*m_data[_j];
  return crossProd;

}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
void
Matrix<T>::print() const
{
  unsigned int i, j;
  T * p;
  for( i=0, p=m_data; i < m_rows; i++ )
  {
    for( j=0; j < m_cols; j++ )
      cout << setw(12) << *p++;
    cout << endl;
  }
}
//-----------------------------------------------------

//-----------------------------------------------------
//template <class T> inline
//const bool
//Matrix<T>::equalWithTolerance( const Matrix<T>& _m, const T& _a ) const
//{
//  assert( !((m_rows!=_m.m_rows) || (m_cols!=_m.m_cols)) );
//  T eps_rel;
//  for( unsigned int z=0; z<m_rows; ++z )
//  {
//    for( unsigned int s=0; s<m_cols; ++s )
//    {
//      eps_rel = Math::relError( m_data[z*m_cols+s], _m.m_data[z*m_cols+s] );
//      // check tolerance
//      if( fabs(eps_rel) > fabs(_a) )
//      {
//        cout << z << ","<< s << ": "<< fabs(eps_rel) << ">" << fabs(_a) << endl;
//        return false;
//      }
//    }
//  }
//  return true;
//}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Vector<T> 
Matrix<T>::solve( Vector<T>& _v )
{
  int i, j, k;
  double m;
  for(j=1; j<m_cols; j++)
  {
    for(i=j; i<m_cols; i++)
    {
      m = m_data[i*m_cols+j-1] / m_data[(j-1)*m_cols+j-1];
      for(k=0; k<m_cols; k++)
      {
        m_data[i*m_cols+k] -= m * m_data[(j-1)*m_cols+k];                
      }
      _v(i) -= m * _v(j-1);
    }
  }
  // Perform back substitution
  Vector<T> x(m_cols);
  x(m_cols-1) = _v(m_cols-1) / m_data[(m_cols-1)*m_cols+(m_cols-1)];
  for( j=m_cols-2; j>=0; --j )
  {
    m = 0;
    for( i=j+1; i<m_cols; i++ )
    {
      m +=  m_data[j*m_cols+i] * x(i);
    }
    x(j) = (_v(j) - m) /  m_data[j*m_cols+j];
  }
  return x;
}
//-----------------------------------------------------



////////////////////////////////////////////////
// GLOBAL OPERATORS
////////////////////////////////////////////////

//-----------------------------------------------------
template <class T> inline
ostream&
operator << (ostream& s, const Matrix<T>& m)
{
  for( unsigned int i=0; i < m.getRows(); i++ )
  {
    for( unsigned int j=0; j < m.getCols(); j++ )
      s << m(i,j) << '\t';
    s << endl;
  }
  return s;
}
//-----------------------------------------------------

//-----------------------------------------------------
template <class T> inline
Matrix<T>
operator*(const T &a, const Matrix<T> &x)
{
  unsigned int i, j;
  unsigned int R = x.getRows();
  unsigned int C = x.getCols();
  Matrix<T> result(R,C);
  for(i=0;i<R;i++)
    for(j=0;j<C;j++)
      result(i,j) = x(i,j)*a;
  return result;
}
//-----------------------------------------------------

//-----------------------------------------------------
// unary transpose operator
template <class T> inline
Matrix<T>
operator ~ (const Matrix<T> &x)
{
  unsigned int i, j;
  unsigned int R = x.getRows();
  unsigned int C = x.getCols();
  Matrix<T> result(C,R);
  for(i=0;i<R;i++)
    for(j=0;j<C;j++)
      result(j,i) = x(i,j);
  return result;
}
//-----------------------------------------------------


#endif // _MATRIX_H_



/*
// unary negation operator
MAT_TEMPLATE inline matrixT
matrixT::operator - () _NO_THROW
{
   matrixT temp(_m->Row,_m->Col);

   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 temp._m->Val[i][j] = - _m->Val[i][j];

   return temp;
}

// binary addition operator
MAT_TEMPLATE inline matrixT
operator + (const matrixT& m1, const matrixT& m2) _THROW_MATRIX_ERROR
{
   matrixT temp = m1;
   temp += m2;
   return temp;
}


// binary power operator
MAT_TEMPLATE inline matrixT
operator ^ (const matrixT& m, const size_t& pow) _THROW_MATRIX_ERROR
{
   matrixT temp = m;
   temp ^= pow;
   return temp;
}

// unary transpose operator
MAT_TEMPLATE inline matrixT
operator ~ (const matrixT& m) _NO_THROW
{
   matrixT temp(m.ColNo(),m.RowNo());

   for (size_t i=0; i < m.RowNo(); i++)
      for (size_t j=0; j < m.ColNo(); j++)
      {
         T x = m(i,j);
	      temp(j,i) = x;
      }
   return temp;
}

// unary inversion operator
MAT_TEMPLATE inline matrixT
operator ! (const matrixT m) _THROW_MATRIX_ERROR
{
   matrixT temp = m;
   return temp.Inv();
}

// inversion function
MAT_TEMPLATE inline matrixT
matrixT::Inv () _THROW_MATRIX_ERROR
{
   size_t i,j,k;
   T a1,a2,*rowptr;

   if (_m->Row != _m->Col)
      REPORT_ERROR( "matrixT::operator!: Inversion of a non-square matrix");

   matrixT temp(_m->Row,_m->Col);
   if (_m->Refcnt > 1) clone();


   temp.Unit();
   for (k=0; k < _m->Row; k++)
   {
      int indx = pivot(k);
      if (indx == -1)
	      REPORT_ERROR( "matrixT::operator!: Inversion of a singular matrix");

      if (indx != 0)
      {
	      rowptr = temp._m->Val[k];
	      temp._m->Val[k] = temp._m->Val[indx];
	      temp._m->Val[indx] = rowptr;
      }
      a1 = _m->Val[k][k];
      for (j=0; j < _m->Row; j++)
      {
	      _m->Val[k][j] /= a1;
	      temp._m->Val[k][j] /= a1;
      }
      for (i=0; i < _m->Row; i++)
	   if (i != k)
	   {
	      a2 = _m->Val[i][k];
	      for (j=0; j < _m->Row; j++)
	      {
	         _m->Val[i][j] -= a2 * _m->Val[k][j];
	         temp._m->Val[i][j] -= a2 * temp._m->Val[k][j];
	      }
	   }
   }
   return temp;
}

// solve simultaneous equation
MAT_TEMPLATE inline matrixT
matrixT::Solve (const matrixT& v) const _THROW_MATRIX_ERROR
{
   size_t i,j,k;
   T a1;

   if (!(_m->Row == _m->Col && _m->Col == v._m->Row))
      REPORT_ERROR( "matrixT::Solve():Inconsistent matrices!");

   matrixT temp(_m->Row,_m->Col+v._m->Col);
   for (i=0; i < _m->Row; i++)
   {
      for (j=0; j < _m->Col; j++)
	 temp._m->Val[i][j] = _m->Val[i][j];
      for (k=0; k < v._m->Col; k++)
	 temp._m->Val[i][_m->Col+k] = v._m->Val[i][k];
   }
   for (k=0; k < _m->Row; k++)
   {
      int indx = temp.pivot(k);
      if (indx == -1)
	 REPORT_ERROR( "matrixT::Solve(): Singular matrix!");

      a1 = temp._m->Val[k][k];
      for (j=k; j < temp._m->Col; j++)
	 temp._m->Val[k][j] /= a1;

      for (i=k+1; i < _m->Row; i++)
      {
	 a1 = temp._m->Val[i][k];
	 for (j=k; j < temp._m->Col; j++)
	   temp._m->Val[i][j] -= a1 * temp._m->Val[k][j];
      }
   }
   matrixT s(v._m->Row,v._m->Col);
   for (k=0; k < v._m->Col; k++)
      for (int m=int(_m->Row)-1; m >= 0; m--)
      {
	 s._m->Val[m][k] = temp._m->Val[m][_m->Col+k];
	 for (j=m+1; j < _m->Col; j++)
	    s._m->Val[m][k] -= temp._m->Val[m][j] * s._m->Val[j][k];
      }
   return s;
}

// set zero to all elements of this matrix
MAT_TEMPLATE inline void
matrixT::Null (const size_t& row, const size_t& col) _NO_THROW
{
    if (row != _m->Row || col != _m->Col)
	realloc( row,col);

    if (_m->Refcnt > 1) 
	clone();

    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = T(0);
    return;
}

// set zero to all elements of this matrix
MAT_TEMPLATE inline void
matrixT::Null() _NO_THROW
{
    if (_m->Refcnt > 1) clone();   
    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
		_m->Val[i][j] = T(0);
    return;
}

// set this matrix to unity
MAT_TEMPLATE inline void
matrixT::Unit (const size_t& row) _NO_THROW
{
    if (row != _m->Row || row != _m->Col)
	realloc( row, row);
	
    if (_m->Refcnt > 1) 
	clone();

    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = i == j ? T(1) : T(0);
    return;
}

// set this matrix to unity
MAT_TEMPLATE inline void
matrixT::Unit () _NO_THROW
{
    if (_m->Refcnt > 1) clone();   
    size_t row = min(_m->Row,_m->Col);
    _m->Row = _m->Col = row;

    for (size_t i=0; i < _m->Row; i++)
	for (size_t j=0; j < _m->Col; j++)
	    _m->Val[i][j] = i == j ? T(1) : T(0);
    return;
}

// private partial pivoting method
MAT_TEMPLATE inline int
matrixT::pivot (size_t row)
{
  int k = int(row);
  double amax,temp;

  amax = -1;
  for (size_t i=row; i < _m->Row; i++)
    if ( (temp = abs( _m->Val[i][row])) > amax && temp != 0.0)
     {
       amax = temp;
       k = i;
     }
  if (_m->Val[k][row] == T(0))
     return -1;
  if (k != int(row))
  {
     T* rowptr = _m->Val[k];
     _m->Val[k] = _m->Val[row];
     _m->Val[row] = rowptr;
     return k;
  }
  return 0;
}

// calculate the determinant of a matrix
MAT_TEMPLATE inline T
matrixT::Det () const _THROW_MATRIX_ERROR
{
   size_t i,j,k;
   T piv,detVal = T(1);

   if (_m->Row != _m->Col)
      REPORT_ERROR( "matrixT::Det(): Determinant a non-square matrix!");
   
   matrixT temp(*this);
   if (temp._m->Refcnt > 1) temp.clone();

   for (k=0; k < _m->Row; k++)
   {
      int indx = temp.pivot(k);
      if (indx == -1)
	 return 0;
      if (indx != 0)
	 detVal = - detVal;
      detVal = detVal * temp._m->Val[k][k];
      for (i=k+1; i < _m->Row; i++)
      {
	 piv = temp._m->Val[i][k] / temp._m->Val[k][k];
	 for (j=k+1; j < _m->Row; j++)
	    temp._m->Val[i][j] -= piv * temp._m->Val[k][j];
      }
   }
   return detVal;
}

// calculate the norm of a matrix
MAT_TEMPLATE inline T
matrixT::Norm () _NO_THROW
{
   T retVal = T(0);

   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 retVal += _m->Val[i][j] * _m->Val[i][j];
   retVal = sqrt( retVal);

   return retVal;
}

// calculate the condition number of a matrix
MAT_TEMPLATE inline T
matrixT::Cond () _NO_THROW
{
   matrixT inv = ! (*this);
   return (Norm() * inv.Norm());
}

// calculate the cofactor of a matrix for a given element
MAT_TEMPLATE inline T
matrixT::Cofact (size_t row, size_t col) _THROW_MATRIX_ERROR
{
   size_t i,i1,j,j1;

   if (_m->Row != _m->Col)
      REPORT_ERROR( "matrixT::Cofact(): Cofactor of a non-square matrix!");

   if (row > _m->Row || col > _m->Col)
      REPORT_ERROR( "matrixT::Cofact(): Index out of range!");

   matrixT temp (_m->Row-1,_m->Col-1);

   for (i=i1=0; i < _m->Row; i++)
   {
      if (i == row)
	continue;
      for (j=j1=0; j < _m->Col; j++)
      {
	 if (j == col)
	    continue;
	 temp._m->Val[i1][j1] = _m->Val[i][j];
	 j1++;
      }
      i1++;
   }
   T  cof = temp.Det();
   if ((row+col)%2 == 1)
      cof = -cof;

   return cof;
}


// calculate adjoin of a matrix
MAT_TEMPLATE inline matrixT
matrixT::Adj () _THROW_MATRIX_ERROR
{
   if (_m->Row != _m->Col)
      REPORT_ERROR( "matrixT::Adj(): Adjoin of a non-square matrix.");

   matrixT temp(_m->Row,_m->Col);

   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	  temp._m->Val[j][i] = Cofact(i,j);
   return temp;
}

// Determine if the matrix is singular
MAT_TEMPLATE inline bool
matrixT::IsSingular () _NO_THROW
{
   if (_m->Row != _m->Col)
      return false;
   return (Det() == T(0));
}

// Determine if the matrix is diagonal
MAT_TEMPLATE inline bool
matrixT::IsDiagonal () _NO_THROW
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=0; i < _m->Row; i++)
     for (size_t j=0; j < _m->Col; j++)
	if (i != j && _m->Val[i][j] != T(0))
	  return false;
   return true;
}

// Determine if the matrix is scalar
MAT_TEMPLATE inline bool
matrixT::IsScalar () _NO_THROW
{
   if (!IsDiagonal())
     return false;
   T v = _m->Val[0][0];
   for (size_t i=1; i < _m->Row; i++)
     if (_m->Val[i][i] != v)
	return false;
   return true;
}

// Determine if the matrix is a unit matrix
MAT_TEMPLATE inline bool
matrixT::IsUnit () _NO_THROW
{
   if (IsScalar() && _m->Val[0][0] == T(1))
     return true;
   return false;
}

// Determine if this is a null matrix
MAT_TEMPLATE inline bool
matrixT::IsNull () _NO_THROW
{
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 if (_m->Val[i][j] != T(0))
	    return false;
   return true;
}

// Determine if the matrix is symmetric
MAT_TEMPLATE inline bool
matrixT::IsSymmetric () _NO_THROW
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 if (_m->Val[i][j] != _m->Val[j][i])
	    return false;
   return true;
}
	   
// Determine if the matrix is skew-symmetric
MAT_TEMPLATE inline bool
matrixT::IsSkewSymmetric () _NO_THROW
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=0; i < _m->Row; i++)
      for (size_t j=0; j < _m->Col; j++)
	 if (_m->Val[i][j] != -_m->Val[j][i])
	    return false;
   return true;
}
   
// Determine if the matrix is upper triangular
MAT_TEMPLATE inline bool
matrixT::IsUpperTriangular () _NO_THROW
{
   if (_m->Row != _m->Col)
      return false;
   for (size_t i=1; i < _m->Row; i++)
      for (size_t j=0; j < i-1; j++)
	 if (_m->Val[i][j] != T(0))
	    return false;
   return true;
}

// Determine if the matrix is lower triangular
MAT_TEMPLATE inline bool
matrixT::IsLowerTriangular () _NO_THROW
{
   if (_m->Row != _m->Col)
      return false;

   for (size_t j=1; j < _m->Col; j++)
      for (size_t i=0; i < j-1; i++)
	 if (_m->Val[i][j] != T(0))
	    return false;

   return true;
}


*/
/*
#define TINY 1.0e-20;
//we assume fabs(ElType) is defined
//assignment of doubles to ElType is defined
template <class ElType>
void ludcmp(matrix<ElType>& a, vector<int>& indx,double& d)
{
	int i,imax,j,k;
	ElType  big,dum,sum,temp;
	int n=a.rows;
	vector<ElType> vv(n);
	assert(a.rows == a.cols);
	d=1.0;
	for (i=0;i<n;i++)
	{
		big=0.0;
		for (j=0;j<n;j++) if ((temp=fabs(a[i][j])) > big) big=temp;
		if (big == 0.0) cerr << "Singular matrix in routine LUDCMP" << endl;
		vv[i]=1.0/big;
	}
	for (j=0;j<n;j++)
	{
		for (i=0;i<j;i++)
		{
			sum=a[i][j];
			for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		big=0.0;
		for (i=j;i<n;i++)
		{
			sum=a[i][j];
			for (k=0;k<j;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*fabs(sum)) >= big)
			{
				big=dum;
				imax=i;
			}
		}
		if (j != imax)
		{
			for (k=0;k<n;k++)
			{
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			d = -(d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]=TINY;
		if (j != n-1) {
			dum=1.0/(a[j][j]);
			for (i=j+1;i<n;i++) a[i][j] *= dum;
		}
	}
}
#undef TINY
template <class ElType>
void lubksb(matrix<ElType>& a,vector<int>& indx,vector<ElType>& b)
{
	int i,ip,j;
	ElType sum;
	int n=a.rows;
	for (i=0;i<n;i++)
	{
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		for (j=0;j<=i-1;j++) sum -= a[i][j]*b[j];
		b[i]=sum;
	}
	for (i=n-1;i>=0;i--)
	{
		sum=b[i];
		for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}

*/
